Lattice density-functional theory on graphene 
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A density-functional approach on the hexagonal graphene lattice is developed using an exact 
numerical solution to the Hubbard model as the reference system. Both nearest-neighbour and up 
to third nearest- neighbour hoppings are considered and exchange-correlation potentials within the 
local density approximation are parameterized for both variants. The method is used to calculate 
the ground-state energy and density of graphene flakes and infinite graphene sheet. The results are 
found to agree with exact diagonalization for small systems, also if local impurities are present. In 
addition, correct ground-state spin is found in the case of large triangular and bowtie flakes out of 
the scope of exact diagonalization methods. 



I. INTRODUCTION 

Density-functional theory (DFT) is one of the most 
widely-used tools in the field of electronic structure sim- 
ulations. Based on the Hohenberg-Kohn theorems^ and 
Kohn-Sham equations,— it allows one to treat interacting 
many-body problems as a collection of one-body prob- 
lems with an effective potential. The many-body effects 
are taken into account through this exchange-correlation 
potential, whose functional form is not known exactly. 

Despite its apparent simplicity, the Hubbard model, 
containing only the on-site contribution of the Coulomb 
interaction, has been found to describe a wide range of 
correlated materials and phenomena.- As the model can- 
not be solved analytically above one dimension, numer- 
ical or approximative methods, such as exact diagonal- 
ization (ED) on small systems, quantum Monte Carlo or 
perturbation theory, have to be applied. For example in 
transport calculations on graphene, the mean-field ap- 
proximation has succesfully been applied to allow treat- 
ment of systems consisting of few hundreds of atomSf 4 - 

The DFT and the Hubbard model can be combined 
in several ways£ The most obvious choices are either to 
determine the model parameters based on first principles 
calculations or to incorporate a Hubbard-type interaction 
into the DFT exchange-correlation functional, resulting 
in the so-called DFT+U method^ The model, however, 
can also be considered as an interesting system on its own 
and studied using DFT in a lattice formulation. This 
approach has previously been chosen by, for instance, 
Capelle et al. (named BA-LDA) ^1 and Gunnarsson and 
Schonhammer (named SOFT)£~— and Schindlmayr and 
Godby^- 

The exchange-correlation (XC) functional, the core of 
the many-body treatment in DFT, is naturally needed 
also in the lattice formulation. In the case of the one- 
dimensional Hubbard model, an exact analytical solution 
based on the Bethe Ansatz is known^ and it can be used 
to parameterize the XC functional.- One-dimensional 
Hubbard chains have been widely studied within this 
framework.— -11 ' 13 The method has been found to accu- 
rately reproduce the exact ground-state energy, also in 
the presence of impurities modelled as on-site energies^ 



Also a related method, density matrix functional the- 
ory, has been applied on the Hubbard model, not only 
for the 1-dimensional chain but also for small clusters 
and square lattice in two and three dimensions J^r— The 
functional in this approach has been based on quantum 
Monte Carlo reference data. In general, the density ma- 
trix functional theory has shown a great potential, be- 
ing applicable even in strongly correlated quantum Hall 
droplets*^ 



To the best of our knowledge, no lattice density- 
functional theory (LDFT) studies have previously been 
performed on the two-dimensional hexagonal lattice, or 
in dimensions above one, in general. The LDFT method 
could be used in transport calculations on graphene in- 
stead of the usual mean-field treatment, thus improving 
the description of correlations for large graphene systems. 
It would also enable one to study dynamic phenomena in 
large systems in the form of time-dependent LDFT. 



Graphene nanoflakcs, or graphene quantum dots have 
been proposed as elements of future nanoelectronics de- 
vices. Especially triangular nanoflakes have been a sub- 
ject of recent research interest. Zigzag-edged flakes have 
been found to show non-trivial spin order and they 
have been proposed to function for instance as logic 
devices in nanoelectronics J£r— Also the effect of edge 
termination on their transport properties has recently 
been studied, spin valve and spin rectification prop- 
erties being reported i2I The transmission through tri- 
angular junctions was found tunable through holes in 
the triangular flakes forming the junction^ They do 
exhibit spin-polarized ground-states, both in density- 
functional theory calculations 1 ^— and in theoretical con- 
siderations within the nearest-neighbour tight-binding 
scheme in the absence of interactions! 23 ' 24 Our pur- 
pose is to apply the Hubbard-based LDFT method 
on large trianglar nanoflakes, also including the up to 
third nearest-neighbour hopping neglected in the earlier 
calculations and show that our method captures the 
essentials of these phenomena. 



II. LATTICE DENSITY-FUNCTIONAL 
THEORY 

In DFT, density is the main variable instead of the 
wave function. In the lattice formulation, the discrete 
density {rii}, also called the site occupations, takes this 
place. The Hohenberg-Kohn theorems apply also in 
the discrete formulation 11 and the Kohn-Sham equations 
read 

HKStpi = £itpi, (1) 

where Hhs, the Kohn-Sham Hamiltonian, is 

Hks = T + V cxt + V H + V xc . (2) 

In these equations, ipi and ti are the Kohn-Sham eigen- 
states and -energies, respectively, T is the kinetic energy 
operator, V ex t the external potential, Vh the Hartree po- 
tential and V xc the exchange-correlation potential. These 
are simply defined on the set of lattice sites, instead of 
as a continuous function of the position variable, in the 
discrete formulation. As the exchange-correlation poten- 
tial is density-dependent, V xc ({rii}), Eqs. (UJ have to be 
solved self-consistently. Inclusion of spin in the model 
leads to separate Kohn-Sham equations for each spin 
species that are coupled through V xc . 
The Hubbard model is given by 

h = X (*«CAj<7 + 44^°-) + ^X ( 3 ) 

where Ci„ (cJ CT ) are the usual annihilation (creation) op- 
erators for spin a, hi a is cj^c^, Uj is the hopping am- 
plitude between sites i and j and U is the strength of 
the on-site interaction between opposite spins. The set 
of hopping amplitudes {Uj} determines the geometry of 
the system. In the case of graphene, hoppings up to the 
first or third neighbour on the hexagonal lattice are rou- 
tinely included, given by either t\ = —2.7 eV in the case 
of nearest-neighbour hopping or t\ = —2.7 cV, t2 = —0.2 
eV and t% — —0.18 eV in the case of up to third- nearest 
neighbour hopping. 4 ' 28 In this work, we choose our unit 
of energy to be t = —t\, scaling the other hopping am- 
plitudes accordingly. 

The kinetic term T of Eq. is given in our lattice for- 
mulation by the hopping term of the Hubbard Hamilto- 
nian. The Hartree potential Vh, of which only the on-site 
contribution is taken into account in the Hubbard model, 
is considered in the sense of the mean-field Hartree-Fock 
approximation (HF), and the corresponding Hamiltonian 
is 

H H f = X Uj{c\ a c.j a +c} a Ci a ) + U^2h i<T n i - (T , (4) 

where rij CT the electron density of electrons with spin a 
on the site i. The second term in Eq. (0| gives the 
interaction contribution. 
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FIG. 1. (a) The 12-atom supercell of the infinite graphene 
sheet, chosen for the determination of the exchange- 
correlation energy and potential. The dashed lines show the 
connections over the periodic boundaries, (b) A Cie flake 
chosen as the first test system of our LDFT method. 

We assume V xc to depend only on the total density 
rii = Titf- + Tin an d n °t separately on the spin-resolved 
densities. In the LDFT calculation, the Kohn-Sham or- 
bitals ipi,r are solved from Eq. ([1]) and the new density is 
calculated from the occupied orbitals of each spin species, 

nia= £ |^ CT | 2 . (5) 

jtrGocc 

The new XC potential is then determined based on the 
occupations and the system is iterated until convergence 
of the site occupations {ni a }. If necessary, the conver- 
gence may be facilitated using a standard mixing proce- 
dure in which the new density is calculated as a linear 
combination of the current and old density. After con- 
vergence, the total energy of the ground-state is obtained 
from the Kohn-Sham eigenenergies and density as 29 

E = X Zir + E ™({ni},U) 

- >, niV xc (ni, U) - U } j njfTHi. (6) 

i i 



III. DETERMINING THE XC FUNCTIONAL 

The exchange-correlation energy and functional con- 
tain the many-body effects of the original problem. In 
the spirit of the XC functional for the ID Hubbard chain 
by Lima et al.^Z- we choose to determine the functional 
within the local density approximation, first introduced 
by Schonhammer et al. in Ref. fiol . The XC energy and 
potential on a given site thus depend only on the local 
density at that site and not, for example, on the density 
difference between the site and its neighbours. The to- 
tal XC energy of the system is then simply obtained by 
summing over the contributions on each of the sites, 

E xc (ni,U) = X e xc(^i, U), (7) 

i 

where e xc (n,i,U) is the XC energy for a homogeneous 
reference system, defined as 
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FIG. 2. (Color online) The XC energy and potential fitted to the infinite graphene sheet with nearest-neighbour hopping for 
S z = 0. Top row: INN hopping, bottom row: 3NN hopping (a/e) and (c/g): numerical XC energy and potential from ED 
(b/f) and (d/h): fitted analytic functions of form Eq. ([9]). 



*.(n,U) 



E(n, U) — E (n, U) 



(8) 



Here, N s is the number of sites in the reference system, 
E(n,U) is the exact reference energy given in our case 
by the exact diagonalization of the Hubbard Hamiltonian 
at filling n = (Nf + Ni)/N s , where N a is the number of 
electrons of the spin species a, and E HF (n,U) is the 
energy of the same system calculated using the Hartrcc- 
Fock Hamiltonian, Eq. ((4]) . For each number of electrons 
+ N±, multiple values for the z-projection of the spin, 
S z = |(iV^ — N±), were considered. More details on the 
exact diagonalization technique can be found in Refs. l30l 
and [M|. 

As the Hubbard model is numerically exactly diago- 
nalizable only up to approximately 16-sites on present 
computers in the general nonsymmetric case, we chose a 
12-atom supercell of the infinite graphene sheet as our 
reference system, shown in Fig. [T] As graphene flakes 
were thought to be the first system to apply the method 
to, no /c-dependence was taken into account in the func- 
tional and the calculations were performed at the T point 
with k = 0, introducing no phase shift to the wave func- 
tion when crossing the supercell boundary. 

The Hubbard model was solved for all combinations 
of Nf and N± for U = — At both by exactly diag- 
onalizing the many-body Hubbard Hamiltonian in the 
Fock basis, and in the single-electron framework using 
the the Hartree-Fock Hamiltonian. The range of the U 
values was chosen to be moderate but to well extend over 
the range » t relevant for graphene calculations. 4 Fig. [5] 
shows the numerical data and the fit for the S z = case 



both for the XC energy and potential, on the top row for 
nearest-neighbour hopping and on the bottom row for up 
to third-nearest neighbour hopping. The XC energy was 
fitted to a function of form 



,(n, U) = ai(e 



1 e 



-(a 3 |n — l|-a 4 ) 2 



(9) 



The values for the coefficients on are given in Table |H 
both in the case of nearest-neighbour and up to third 
nearest-neighbour hopping. The XC potential was then 
obtained as 



<9e x 



On 



(10) 



To be rigorous, a functional derivative should be taken 
from e xc . In this case, however, this reduces to the usual 
derivative. In general the fit is very good, although we 
note that the largest deviations occur for U > 3t near and 
at half- filling, n = 1. We also note that the derivative 
discontinuity of the XC energy, seen also in Fig. [2j;dgh, 
is correctly contained in our functional. 

The inclusion of the up to third-nearest neighbour 
(3NN) hopping changes the form of the XC energy sur- 
face only slightly (Fig. [21 bottom row). The particle- 
hole symmetry of the nearest-neighbour (INN) Hubbard 
model is lost through the inclusion of the further hop- 
pings. This causes the XC energy surface to be no longer 
symmetric about n = 1, leading to a separate set of coef- 
ficients ctfj for n < 1 and n > 1. The coefficients for n > 1 
are given in TableUas a*. Comparing the coefficients for 
the INN and 3NN cases, we see that the magnitude of 
the XC energy is slightly smaller for 3NN hopping as the 
coefficient a\ mainly determines this magnitude. Also, 
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TABLE I. The coefficients a % for Eq. ©. INN stands for 
nearest-neighbour hopping and 3NN for up to third nearest- 
neighbour hopping. In the 3NN case the coefficients are given 
separately for n < 1 (a,) and n > 1 (a*). Coefficients based 
on both S z = and Sz = \ reference data are given. See text 
for details. 

INN INN 3NN 3NN 

S z S z 2 S z — S z — 2 



O'l 

„ * 
Oil 

OL2 

Q*2 

0L3 

«3 

CX4 

at 



0.22554 0.24732 



0.05776 0.04919 



2.67381 2.70556 



2.02445 2.22949 



0.17094 
0.17458 
0.08322 
0.08657 
2.77940 
2.73657 
2.09609 
2.02992 



0.15972 
0.15071 
0.09122 
0.10440 
2.98123 
2.82716 
2.26728 
2.10459 



for 3NN hopping the values of the coefficients for n < 1 
and n > 1 are very close to each other. 

The XC energy surfaces for higher values of S z — 
§>•■• were also calculated. As the number of ac- 
cessible density values and thus reference data points de- 
creases with increasing S z due to the finite supercell, the 
fitting of the XC energy surfaces become more and more 
ambiguous. In addition, in Section IIVI we demonstrate 
that our approach captures well locally spin-polarized 
systems using only the S z = potential, although the 
parameters for the S z = \ potential were also determined 
and are shown in Table [ 



IV. RESULTS 

A. Comparing LDFT to exact diagonalization 

As a first test for our functional, we compare the 
ground-state energy of a 16-atom graphene flake (Fig. Q} 
calculated by exactly diagonalizing the Hubbard Hamil- 
tonian, from the self-consistent solution of the Hartree- 
Fock Hamiltonian, Eq. (j4]), andusing our LDFT method. 
We calculate the system at three different fillings: at 
quarter-filling (Nf = N± = 4) , slightly off quarter-filling 
(JV t = 6, Ni = 5) and at half-filling (JV t = JVj. = 8), 
including only the nearest-neighbour hopping. Two dif- 
ferent S z states, S z — and S z = \ are thus included in 
this comparison. Additionally, we perturb the system by 
adding an on-site energy of magnitude e = — \t\ on the 
two middle sites of the structure. In general, introduc- 
tion of the on-site energies gives to Eqs. ^ and (0} a 
term of the form 



^2 e kt\ a t k<J , 

/cer£ {on — site} 



(ii) 



where the sum runs over the site with an on-site energy 
and efc gives the magnitude of this energy. 
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(d) 



FIG. 3. (Color online) A comparison between the ED, HF 
and LDFT ground-state energies of the Ci6 flake (Fig. {TJ. 
(a) At quarter-filling, (b) With S z = ± (N t = 6, JVj. = 5). 
LDFT calculation performed both using a potential fitted for 
S z = and S z = o reference systems, (c) At quarter-filling 
with an added on-site impurity potential e = —\t\ on the two 
middle sites of the structure (d) At half-filling. The insets 
show the relative errors (in percentage) with respect to the 
ED energy. 
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Fig. [3] shows the ground-state energy as a function of 
U calculated using the three methods. The agreement 
between the LDFT and ED, apart from the half- filled 
case, is very good. This is also illustrated in the insets, 
which show the relative error in percentage with respect 
to the ED energy that is below 1 % away from half-filling. 
At half- filling the agreement is worst, the relative error 
being * 6 %, at U ~ it. This disagreement is partly 
due to the worse fit of the XC energy at half-filling. It 
is also worth noting that the difference between the two 
parameterizations of the potential for S z — and S z = i 
systems is actually in favour of using the S z = potential 
for the S z — | test system (see Fig. [3}d) but this is likely 
to be due to a coincidental cancellation of errors. On the 
other hand, this agreement is also to be expected as the 
V xc fits are very close to each other for these two spin 
states, see Table UJ 

A comparison of the ground-state densities further re- 
solves the issue of poorer results at half-filling. For the 
system without any on-site impurities, the exact den- 
sity is uniform also for non-zero U. The Hartree-Fock 
solution, on the other hand, has an antiferromagnetic 
ground-state density with local spin polarizations. This 
characteristic of the Hartree-Fock solution is also con- 
tained in the LDFT density. In the half-filled case at 
U = 3t, the spin-resolved LDFT density ranges from 
0.20 to 0.80 with the maximal local spin polarization 
5*' max = 0.3 for both spin up and spin down species. 
This separation of spin up and spin down densities is 
also seen in the quarter-filled case but with a smaller 
amplitude due to the lower average density. Thus, al- 
though the LDFT method corrects the ground-state en- 
ergy quite accurately, the method fails to reproduce the 
homogenous ED density. This problem is due to the 
single-configuration wave functions in DFTi^ 

The effect of local spin polarizations can be studied by 
introducing a spin-dependent on-site energy as a pertur- 
bation to the model. This causes the ground-state den- 
sity of the ED solution to become non-uniform. In the 
case of a quarter-filled infinite sheet with the 12-atom 
supercell (Fig. Q}, Nf = N± = 3, calculated including up 
to third nearest-neighbour hoppings with an on-site en- 
ergy eo = — 2\t\ applied on a single site for only the spin 
up species, the ED ground state exhibits spin-dependent 
occupations ranging from 0.14 to 0.60 at U = and from 
0.13 to 0.63 at U = St. As also the ED density is stag- 
gered, the LDFT density agrees better with it, the LDFT 
density varying at U = 3t from 0.11 to 0.65. Fig. 0^ 
shows the energy as a function of U at two values of the 
on-site energy, and Fig. |4Jd the ground-state energy as a 
function of the on-site perturbation strength at two dif- 
ferent non-zero values of U. The LDFT energies agree 
again very well with the exact solutions. 

To conclude this discussion, we note that the compu- 
tational effort of LDFT is roughly the same as HF. In 
comparison to ED, LDFT scales polynomially instead of 
exponentially with the system size. 




(a) 




(b) 



FIG. 4. (Color online) The ground-state energy of an in- 
finite graphene sheet (12-atom supercell) at quarter-filling, 
perturbed with a single spin-dependent on-site energy, using 
ED (full line), HF (dashed line) and third-nearest-neighbour 
LDFT (dash-dotted line). In (a) the strength of the on-site 
energy eo is varied at fixed U and in (b) the on-site energy is 
fixed and U varied. 



B. Triangular flakes 

Proceeding away from the simple, exactly solvable sys- 
tems, we apply the LDFT method on triangular graphene 
nanoflakes. First we study a still small triangle consist- 
ing of 22 sites that has S z = 1 in its ground-state at 
half-filling j 23 ' 24 i 33 As a reference, we compare our LDFT 
calculation to a method we call "partial diagonalization" 
(PD) in which instead of diagonalizing the many-body 
Hamiltonian Eq. ([3]) in the full Hilbert space, we only 
use the space spanned by the non-interacting ground- 
state and its up to double excitations. This method is 
thus equivalent to the singles-doubles configuration inter- 
action method used in quantum chemistry. We calculate 
the non-interacting eigenstates of the system and occupy 
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them according to the filling of the flake. We then singly 
and doubly excite the the system and calculate the ele- 
ments of the Hamiltonian in the set of many-body states 
formed by the non-interacting ground-state and its exci- 
tations. This method is accurate at low values of U up to 
U ~ 2t, ensured by comparing to ED calculations for 6- 
and 10-atom flakes. Flake sizes up to approximately 40 
sites are easily accessible within this method. It should, 
however, be noted that the accuracy of PD at a fixed 
number of included excitations is expected to be reduced 
for larger systems as the correlation energy does not grow 
linearly with the system size. 

Fig. [5] compares the ground-state energy of the dif- 
ferent S z states of the triangular C22 flake calculated 
using PD and nearest-neighbour LDFT using only the 
S z = potential. In Fig. [5^, we see that despite the 
use of the S z = potential, LDFT yields the correct S z 
value for the lowest-energy state at all values of U. The 
inset shows the deviation in percentage from the S z = 
state calculated using PD. At U > 3i the PD energy be- 
comes nearly linear in U and we expect the LDFT energy 
to better describe the physics of the system, as for very 
large U the ground-state should exclude doubly-occupied 
sites and the energy should saturated This linearity in 
the PD solution at larger U leads also to apparently large 
deviations in the inset of Fig. [5^. Figs. [5p and [5J: show 
the ground-state density from PD and LDFT, respec- 
tively, for the lowest-energy state (S z = 1) at U = t and 
Figs. [SJl and [SJ3 at U — 2.8t. The excess spin is local- 
ized on one of the sublattices and mostly on the outmost 
zigzag sites, a feature correctly captured by the LDFT 
solution and also reported in DFT calculations J^— For 
U — t, the deviation from the PD density is hardly no- 
ticeable, whereas for U — 2.8t the deviation is clearly 
seen. Even at strong interaction, however, the qualitative 
features of sublattice-dependent polarization and spin lo- 
calization on the outmost sites still remain. 

In the 3NN triangular flakes, unlike the INN 
case ; 23 ! 24 ' 33 there are no general predictions for the 
ground-state spin. As the hopping within the sublat- 
tice is relatively weak (t% ~ 0.07i), it is not likely that 
the situation changes drastically from the system with 
only nearest-neighbour hopping. Third-nearest neigh- 
bour hopping again connects the two sublattices. Ac- 
cording to a theorem by Lieb^ valid for Hubbard mod- 
els on bipartite lattice with a repulsive interaction at 
half-filling, a 286-atom zigzag-edged triangle having 15 
hexagons in the base row should have S z = 7 in its 
ground-state. The spin is given by the sublattice inbal- 
ance, S z = MNa— Nb\, where N^b) signifies the number 
of sites on sublattice A(B). This is due to the fact that 
in the absence of interaction, the highest occupied single- 
electron state is 14-fold degenerate. Our LDFT approach 
indeed finds the state S z = 7 to be of lowest energy for 
U > 0M. Due to the degeneracy, however, convergence 
problems are encountered at U < OAt as the numerical 
solution is very unstable for the open shell iteration for 
lower values of S z , and the lower-spin states appear to 
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FIG. 5. (Color online) Comparison of PD and LDFT within 
the nearest-neighbour hopping approach for the C22 flake hav- 
ing a spin-polarized ground-state at half-filling. In (a) the 
comparison of the ground-state energy for different values of 
S x . Solid line: LDFT, dotted line: PD. The inset shows the 
deviation in percentage from the PD S z — energy. In (b) 
and (c), the ground-state density at U = t from PD and 
LDFT, respectively. In (d) and (e), the density at U — 2.8t. 
Red (light grey) and blue (dark) signify the up and down spin 
densities, respectively. Half-filling corresponds to the radius 
of the blue circles in the lower left corner in (b) and (d) and 
to the single black circle in (c) and (d), drawn to facilitate 
comparison. 



be energetically more favourable. 

In Fig.[6l we show the energy gap at U = t between the 
lowest-energy state of the 286-atom flake and the other 
S z values at and slightly below half-filling, the total num- 
ber of electrons ranging from 273 to 286. The value of U 
was chosen to be relevant to graphene calculations! 4 * 2 ^ At 
half-filling, the S z — 7 has the lowest energy as expected 
from the Lieb theorem. Slightly off half-filling the ener- 
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FIG. 6. The energy gap between the lowest-energy state and 
the other S z states for the 286-atom triangular graphene flake 
at U = t with different total numbers of electrons, N c \ = 
N-f + N±. As the gap jumps when crossing the diagonal, the 
gap is shown only up to the S z value with minimal energy. 
The x-axis shows the difference in S z from the minimal spin, 
S z>m in = for even N ei and S z ,min = 1/2 for uneven N ci . 



getically favourable S z decreases linearly with the total 
number of electrons but at the same time the magnitude 
of the gap decreases. Even though the states seem to lie 
very close to each other in energy, the order of magnitude 
of the gap at half-filling, O.lt « 0.27eV, transformed into 
temperature, rj 3100K, is high enough for us to predict 
that the high-spin ground-state should be observable in 
room temperature. Away from half-filling, on the other 
hand, the gap magnitude in temperature is of the or- 
der of room temperature, making the system an unlikely 
candidate for practical spintronics applications. 

An estimate for the accuracy of LDFT in determining 
the gap magnitude is obtained by considering the gap in 
C 22 (Fig. P. The gap between the S z = 1 and S z = 
ground-states is 0.13 eV in PD and 0.06 eV in LDFT. 
LDFT thus seems to somewhat underestimate the gap 
but this only supports the conclusion on the rigidity of 
the ground-state against thermal fluctuations. As the 
C22 is the smallest flake that has non-minimal spin in 
its ground-state, the accuracy of the LDFT predictions 
in these systems can not be compared to ED. To com- 
pare the gap in a diagonalizable system, the gap between 
the two lowest states (S z — 1/2 and S z = 3/2) in the 
smallest triangular flake, C13, at half-filling was calcu- 
lated using ED, PD, and LDFT at U = t including up 
to third-nearest neighbour hopping. The gap magnitude 
in ED, 1.52 t, is thus an order of magnitude larger than 
in the spin-polarized flakes, and both PED and LDFT 
overestimate it but only by 3.3% and 4.2%, respectively. 

In order to show that the LDFT method is able to pre- 
dict the spin densities also for more complex geometries, 
Fig. [7] shows the ground-state density of a bowtie-shaped 
flake obtained by combining two triangular flakes. Hop- 
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FIG. 7. (Color online) The ground-state density (S z = 0) 
of a 78-atom bowtie flake at half-filling, previously studied by 
Wang et al. in Ref.lla. Up to third nearest-neighbour hopping 
is included and the value of U is t. Red (light grey) and 
blue (dark) correspond to spin up and spin down densities, 
respectively. 



ping up to third-nearest neighbours is included and the 
flake is calculated at half-filling and U = t. Unlike in 
the triangular flakes, the sublattice imbalance is zero in 
these structures and thus the ground-state has S z = 0. 
Nevertheless, the spin structure is nontrivial as the spin 
up and spin down densities concentrate on the opposite 
sides of the structure. These structures have been previ- 
ously proposed to function as spin logic devices by Wang 
et a/.jis and a comparison with their density- functional 
calculation shows that LDFT indeed captures this non- 
trivial spin order. 

Based on the results of this section, we conclude that 
the LDFT approach captures also the essential features 
of systems with a non-zero total spin in the ground-state 
or with large local spin polarizations. 



V. CONCLUSIONS 

Lattice-density functional theory is a promising candi- 
date to be used instead of the mean-field Hartree-Fock 
approximation for graphene systems too large for ex- 
act methods. The exact ground-state energies are accu- 
rately reproduced, and correct spin properties for large 
graphene flakes are found. Local density approximation 
(LDA) is enough to capture these spin properties and 
no need for an explicitly spin-dependent potential was 
found. The ground-state densities are antiferromagnetic 
with local spin-polarization like the HF solution and the 
computational effort is of the same order as in the mean- 
field Hartree-Fock approximation. 

The method could be applied to transport calculations, 
for instance for graphene nanoribbons. Also, the method 
could be extended into time domain in the form of an 
adiabatic time-dependent lattice density-functional the- 
ory and used to study dynamic phenomena such as the 
possibility of spin-charge separation above one dimen- 
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